\[ \hat{\beta}=(X'X)^{-1}X'Y \] - Se tiene un problema en cuanto a la transpuesta de la matriz \((X'X)\)
Causas - ¿Cuáles son sus consecuencias prácticas?
Incidencia en los errores estándar y sensibilidad
Pruebas
donde
\[ Y_{i} = \beta_{0}+\beta_{1}X_{1}+\beta_{2}X_{2}+\beta_{3}X_{3} + \beta_{4}X_{4}+\beta_{5}X_{5}+u_{i}\nonumber \]
Agreguemos el tiempo: - Las correlaciones muy altas también suelen ser síntoma de multicolinealidad
ajuste.2 <- lm(Y~X1+X2+X3+X4+X5+TIME)
summary(ajuste.2)
##
## Call:
## lm(formula = Y ~ X1 + X2 + X3 + X4 + X5 + TIME)
##
## Residuals:
## Min 1Q Median 3Q Max
## -381.7 -167.6 13.7 105.5 488.9
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.727e+04 2.324e+04 2.895 0.02005 *
## X1 -2.051e+00 8.710e+00 -0.235 0.81974
## X2 -2.733e-02 3.317e-02 -0.824 0.43385
## X3 -1.952e+00 4.767e-01 -4.095 0.00346 **
## X4 -9.582e-01 2.162e-01 -4.432 0.00219 **
## X5 5.134e-02 2.340e-01 0.219 0.83181
## TIME 1.585e+03 4.827e+02 3.284 0.01112 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 295.6 on 8 degrees of freedom
## Multiple R-squared: 0.9955, Adjusted R-squared: 0.9921
## F-statistic: 295.8 on 6 and 8 DF, p-value: 6.041e-09
cor(cbind(X1,X2,X3,X4,X5,TIME))
## X1 X2 X3 X4 X5 TIME
## X1 1.0000000 0.9936689 0.5917342 0.4689737 0.9833160 0.9908435
## X2 0.9936689 1.0000000 0.5752804 0.4587780 0.9896976 0.9947890
## X3 0.5917342 0.5752804 1.0000000 -0.2032852 0.6747642 0.6465669
## X4 0.4689737 0.4587780 -0.2032852 1.0000000 0.3712428 0.4222098
## X5 0.9833160 0.9896976 0.6747642 0.3712428 1.0000000 0.9957420
## TIME 0.9908435 0.9947890 0.6465669 0.4222098 0.9957420 1.0000000
Regresemos una de las variables
ajuste.3<- lm(X1~X2+X3+X4+X5+TIME)
summary(ajuste.3)
##
## Call:
## lm(formula = X1 ~ X2 + X3 + X4 + X5 + TIME)
##
## Residuals:
## Min 1Q Median 3Q Max
## -18.8602 -4.3277 -0.3175 4.3726 14.8438
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.529e+03 7.288e+02 2.098 0.0653 .
## X2 2.543e-03 9.453e-04 2.690 0.0248 *
## X3 3.056e-02 1.514e-02 2.019 0.0742 .
## X4 1.011e-02 7.559e-03 1.337 0.2140
## X5 -1.263e-02 7.903e-03 -1.598 0.1445
## TIME -1.621e+01 1.766e+01 -0.918 0.3826
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11.31 on 9 degrees of freedom
## Multiple R-squared: 0.9923, Adjusted R-squared: 0.9881
## F-statistic: 232.5 on 5 and 9 DF, p-value: 3.127e-09
tolerancia<-1-0.9923
Factor de inflación de la varianza
Si este valor es mucho mayor que 10 y se podría concluir que si hay multicolinealidad
vif <- 1/tolerancia
vif
## [1] 129.8701
Ahora vamos a usar el paquete AER:
library(AER)
vif1 <- vif(ajuste.2)
Raux <- (vif1-1)/vif1
Rglobal <- 0.9955
Rglobal-Raux
## X1 X2 X3 X4 X5
## 0.003181137 -0.003829181 0.026533869 0.254649059 -0.001623122
## TIME
## -0.003160352
Se podría no hacer nada ante este problema. O se puede tratar con transformaciones. Deflactamos el PIB: PIB_REAL <- X2/X1
# La variable X5 (población)
# esta correlacionada con el tiempo
PIB_REAL <- X2/X1
ajuste.4<-lm(Y~PIB_REAL+X3+X4)
summary(ajuste.4)
##
## Call:
## lm(formula = Y ~ PIB_REAL + X3 + X4)
##
## Residuals:
## Min 1Q Median 3Q Max
## -760.29 -197.71 -53.69 234.77 603.15
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 42716.5646 710.1206 60.154 3.31e-15 ***
## PIB_REAL 72.0074 3.3286 21.633 2.30e-10 ***
## X3 -0.6810 0.1693 -4.023 0.00201 **
## X4 -0.8392 0.2206 -3.805 0.00292 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 389 on 11 degrees of freedom
## Multiple R-squared: 0.9893, Adjusted R-squared: 0.9864
## F-statistic: 339.5 on 3 and 11 DF, p-value: 4.045e-11
vif(ajuste.4)
## PIB_REAL X3 X4
## 3.054580 2.346489 2.318500
ajuste.5<-lm(Y~PIB_REAL+X3+X4)
summary(ajuste.5)
##
## Call:
## lm(formula = Y ~ PIB_REAL + X3 + X4)
##
## Residuals:
## Min 1Q Median 3Q Max
## -760.29 -197.71 -53.69 234.77 603.15
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 42716.5646 710.1206 60.154 3.31e-15 ***
## PIB_REAL 72.0074 3.3286 21.633 2.30e-10 ***
## X3 -0.6810 0.1693 -4.023 0.00201 **
## X4 -0.8392 0.2206 -3.805 0.00292 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 389 on 11 degrees of freedom
## Multiple R-squared: 0.9893, Adjusted R-squared: 0.9864
## F-statistic: 339.5 on 3 and 11 DF, p-value: 4.045e-11
vif(ajuste.5)
## PIB_REAL X3 X4
## 3.054580 2.346489 2.318500
Ocurre cuando la varianza no es constante.
¿Cuál es la naturaleza de la heterocedasticidad?
Método gráfico
Veamos las pruebas de detección en un ejemplo
Abrir la base de datos wage1 de Wooldrigde
Hacer un gráfico de los valores estimados y los residuos al cuadrado
bptest(objeto): si el pvalor es inferior a \(0.05\), Ho: HomocedasticidadEl códgio en R sería:
ajuste1 <- lm(lwage~casados+casadas+solteras+educ+exper+
expersq+tenure+tenursq)
summary(ajuste1)
##
## Call:
## lm(formula = lwage ~ casados + casadas + solteras + educ + exper +
## expersq + tenure + tenursq)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.89697 -0.24060 -0.02689 0.23144 1.09197
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.3213780 0.1000090 3.213 0.001393 **
## casados 0.2126756 0.0553572 3.842 0.000137 ***
## casadas -0.1982677 0.0578355 -3.428 0.000656 ***
## solteras -0.1103502 0.0557421 -1.980 0.048272 *
## educ 0.0789103 0.0066945 11.787 < 2e-16 ***
## exper 0.0268006 0.0052428 5.112 4.50e-07 ***
## expersq -0.0005352 0.0001104 -4.847 1.66e-06 ***
## tenure 0.0290875 0.0067620 4.302 2.03e-05 ***
## tenursq -0.0005331 0.0002312 -2.306 0.021531 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3933 on 517 degrees of freedom
## Multiple R-squared: 0.4609, Adjusted R-squared: 0.4525
## F-statistic: 55.25 on 8 and 517 DF, p-value: < 2.2e-16
residuos <- resid(ajuste1)
sqresid <- residuos^2
y_techo <- fitted(ajuste1)
plot(y_techo,sqresid)
plot(fitted(ajuste1),resid(ajuste1))
# Usando el "default" de R:
par(mfrow=c(2,2))
plot(ajuste1)
par(mfrow=c(1,1))
library(sandwich)
library(lmtest)
#install.packages("lmSupport")
library(lmSupport)
# Test para ver si hay heterocedasticidad
residuos <- resid(ajuste1)
sqresid <- (residuos)^2
ajuste2 <- lm(sqresid~casados+casadas+solteras+educ+exper+expersq+tenure+tenursq)
summary(ajuste2)
##
## Call:
## lm(formula = sqresid ~ casados + casadas + solteras + educ +
## exper + expersq + tenure + tenursq)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.2346 -0.1237 -0.0887 0.0202 3.4689
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.029e-02 6.893e-02 0.729 0.46603
## casados -4.870e-02 3.816e-02 -1.276 0.20241
## casadas -5.147e-02 3.986e-02 -1.291 0.19727
## solteras 4.162e-03 3.842e-02 0.108 0.91379
## educ 3.849e-03 4.614e-03 0.834 0.40462
## exper 1.008e-02 3.614e-03 2.790 0.00546 **
## expersq -2.071e-04 7.611e-05 -2.720 0.00674 **
## tenure 4.763e-04 4.661e-03 0.102 0.91864
## tenursq 8.670e-05 1.594e-04 0.544 0.58672
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2711 on 517 degrees of freedom
## Multiple R-squared: 0.02507, Adjusted R-squared: 0.009989
## F-statistic: 1.662 on 8 and 517 DF, p-value: 0.105
# F =1.662 y pvalue=0.105 NO EXISTE HETEROCEDASTICIDAD
#Breusch-Pagan test
'bptest es igual a hettest en STATA'
## [1] "bptest es igual a hettest en STATA"
bptest(ajuste1)
##
## studentized Breusch-Pagan test
##
## data: ajuste1
## BP = 13.189, df = 8, p-value = 0.1055
Autocorrelación: correlación entre miembros de series de observaciones ordenadas en el tiempo [como en datos de series de tiempo] o en el espacio [como en datos de corte transversal]:
\[ E(u_i,u_j) \neq 0 \\ i \neq j \] El supuesto es: \[ cov(u_i,u_j|x_i,x_j) = E(u_i,u_j) = 0 \\ i \neq j \]
Método gráfico
Veamos las pruebas de detección en un ejemplo
Abrir la tabla 12.4. Veamos los datos en forma gráfica, y corramos el modelo:
#Indice de compensacion real (salario real)
plot(X,Y)
ajuste.indice<-lm(Y~X)
summary(ajuste.indice)
##
## Call:
## lm(formula = Y ~ X)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.138 -2.130 0.364 2.201 3.632
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 29.5192 1.9424 15.20 <2e-16 ***
## X 0.7137 0.0241 29.61 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.676 on 38 degrees of freedom
## Multiple R-squared: 0.9584, Adjusted R-squared: 0.9574
## F-statistic: 876.5 on 1 and 38 DF, p-value: < 2.2e-16
Revisemos si hay autocorelación:
residuos<- resid(ajuste.indice)
plot(residuos,t="l",xlab="Tiempo")
par(mfrow = c(2,2))
plot(ajuste.indice)
par(mfrow = c(1,1))
Veamos si se trata de una función cuadrática y cúbica
ajuste2 <- lm(Y~X+I(X^2))
summary(ajuste2)
##
## Call:
## lm(formula = Y ~ X + I(X^2))
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.58580 -0.76248 0.09209 0.68442 2.63570
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.622e+01 2.955e+00 -5.489 3.09e-06 ***
## X 1.949e+00 7.799e-02 24.987 < 2e-16 ***
## I(X^2) -7.917e-03 4.968e-04 -15.936 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9669 on 37 degrees of freedom
## Multiple R-squared: 0.9947, Adjusted R-squared: 0.9944
## F-statistic: 3483 on 2 and 37 DF, p-value: < 2.2e-16
ajuste3 <- lm(Y~X+I(X^2)+I(X^3))
summary(ajuste3)
##
## Call:
## lm(formula = Y ~ X + I(X^2) + I(X^3))
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.63265 -0.79419 0.06568 0.66627 2.43810
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.222e+01 1.344e+01 -1.653 0.107060
## X 2.196e+00 5.466e-01 4.018 0.000286 ***
## I(X^2) -1.119e-02 7.178e-03 -1.559 0.127658
## I(X^3) 1.398e-05 3.054e-05 0.458 0.649958
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9774 on 36 degrees of freedom
## Multiple R-squared: 0.9947, Adjusted R-squared: 0.9943
## F-statistic: 2272 on 3 and 36 DF, p-value: < 2.2e-16
Nos quedamos con el ajuste 2
El gráfico de los val ajustados, muestra que se ha eliminado el patron inicial
par(mfrow = c(2,2))
plot(ajuste2)
par(mfrow = c(1,1))
residuos2 <- resid(ajuste2)
plot(residuos2,t="l",xlab="Tiempo")
points(residuos2)
abline(h=0,col="blue")
Cómo debe ser el gráfico
aleatorios=rnorm(40,0,1)
plot(aleatorios,t="l",xlab="Tiempo")
points(aleatorios)
abline(h=0,col="blue")
¿Se parece?
Ejemplo: Pruebas
Ho: No hay autocorrelación
dwtest(ajuste2)
##
## Durbin-Watson test
##
## data: ajuste2
## DW = 1.03, p-value = 0.0001178
## alternative hypothesis: true autocorrelation is greater than 0
¿Cuál es la conclusión?
Otra prueba:
# Ajuste Breuch Godfrey (Ho: No hay autocorrelación)
bgtest(ajuste2,order=4)
##
## Breusch-Godfrey test for serial correlation of order up to 4
##
## data: ajuste2
## LM test = 14.945, df = 4, p-value = 0.004817